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Abstract 

In this paper, we study the asyniptotics and fast computation of the one-sided 
oscillatory Hilbert transforms of the form 

ff+(/(t)e'"*)(x) = / e"^'^^dt, w > 0, a; > 0, 

where the bar indicates the Cauchy principal value and / is a real- valued function 
with analytic continuation in the first quadrant, except possibly a branch point of 
algebraic type at the origin. When x — 0, the integral is interpreted as a Hadamard 
finite-part integral, provided it is divergent. Asymptotic expansions in inverse pow- 
ers of w are derived for each fixed x > 0, which clarify the large w behavior of this 
transform. We then present efficient and affordable approaches for numerical evalu- 
ation of such oscillatory transforms. Depending on the position of x, we classify our 
discussion into three regimes, namely, x = C(l)ora;3>l,0<a:^l and a; = 0. 
Numerical experiments show that the convergence of the proposed methods greatly 
improve when the frequency lo increases. Some extensions to oscillatory Hilbert 
transforms with Bessel oscillators are briefly discussed as well. 
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1 Introduction 



Finite Fourier integrals of the form 



(1.1) 



J a 



with oj > and f{x), g{x) being sufficiently smooth functions have long been the subject 
of intensive study due to their frequent occurrences in wide fields ranging from quan- 
tum chemistry, image analysis, electrodynamics and computerized tomography to fluid 
mechanics |15] . One difficulty in computing integrals (jl.l|) is that, for large frequency 
w, the classical integration methods like Gauss quadrature are inapplicable, since they 
often require many function evaluations which make them highly time consuming. To 
overcome this difficulty, many efficient approaches have been developed and significant 
progress has occurred over the past few years. For instance, based on asymptotic expan- 
sions of (jl.ip as CO tends to infinity, Iserles and N0rsett [151 ES] proposed the asymptotic 
and Filon-type methods to evaluate oscillatory integrals numerically. Other efficient ap- 
proaches include Levin-type methods, numerical steepest descent methods, generalized 
quadrature rules, GMRES methods, modified Clenshaw-Curtis methods, etc.; we refer 
tolalElElEllIlEOlEIllSlllSOlEHlEHland references therein for more information. All 
these methods complement each other but share the advantageous property that their 
accuracy improves greatly when cj increases. 

Recently, oscillatory Hilbert transforms of the form 



have received considerable attention as well. Here, F is an oriented curve in the complex 
plane, / is a complex-valued function satisfying a Holder condition, and the bar denotes 
Cauchy principal value. The interest in theoretical and numerical study of such integral 
transforms arises from the fact that many problems encountered in practice can be repre- 
sented by an integral equation with an oscillatory kernel having a singularity of Cauchy 
type [5l[10l[18]; see also [T7] for numerous applications of Hilbert transforms in applied 
sciences. Although the oscillatory Hilbert transforms (|1.2p bear some resemblances with 
(jl.ip . nevertheless, the singularity of Cauchy type suggests special treatments. It will be 
especially interesting to see, as pointed out in [2S], if the aforementioned methods can 
be extended to handle oscillatory Hilbert transforms. 

For F = [—1, 1], we obtain from (|1.2p the finite oscillatory Hilbert transforms: 



with —l<x<l. An asymptotic expansion of ()1.3p in inverse powers of uj was estab- 
lished by Lyness in [22] based on analytic continuation. Meanwhile, there are several 
numerical schemes available to calculate (jl.Sp , most of which are typically based on in- 
terpolatory type techniques. For example, Okecha [26] proposed to compute (jl.Sp by 




(1.2) 




(1.3) 
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using the Lagrange interpolation polynomial of degree n interpolating / at the n + 1 
zeros of the Legendre polynomial. Capobianco and Criscuolo introduced a numerically 
stable procedure in [5], which relies on an interpolatory procedure at the zeros of Jacobi 
polynomials. In a recent paper [33], Wang and Xiang have presented an integration rule 
of interpolatory type with the aid of the Chebyshev points of the second kind. The rule 
is uniformly convergent with respect to the pole x when / is analytic in a neighborhood 
of the interval [—1, 1], and it can be implemented by means of the fast Fourier transform 
(FFT). If the function / is analytic in a sufficiently large region of the complex plane 
containing [—1,1], then the complex integration method [2^ and the numerical steep- 
est descent method [13] can be extended to compute such integrals efficiently, and the 
accuracy improves greatly as uj increases; see [S^ for details. 

In this paper, we are concerned with one-sided oscillatory Hilbert transforms on the 
positive real axis: 

H+{f{t)e''^'){x) := /°°e*^*^^(ii = lim ( T \ H ^ e'""' ^^dt, x > 0, (1.4) 



t-x e-^0+ \Jn J^+J t 



i.e., r = M"*" in ()1.2p . Here, / is a real- valued function satisfying some conditions. In 
particular, it has an analytic continuation in the first quadrant of the complex plane, 
except possibly a branch point of algebraic type at the origin. When x = 0, the integral 
is interpreted as a Hadamard finite-part integral, provided it is divergent; see Section 
2.21 below for a definition. We point out that the integral ()1.4p is also closely related to 
infinite oscillatory Hilbert transforms on the real axis given by 

'iojt\/^\ ._ f Auit /(O j+ _ i- if , f \ „ioJt /(^) 



Hif{t)e"^'){x) := f e'^'^^^dt = lim / + / e*^*^^dt, x G R. 

Indeed, by assuming x > 0, it is easily seen that 

H{f{t)e'^%x)= r -e-'^'^-^^dt + H+{f{t)e:^%x). (1.5) 
JO t + X 

The first integral on the right hand side of (jl.Sp is the Stieltjes transform of —e~'^^^f{—t), 
which is a regular integral for x > 0. 

A large x expansion of one-sided oscillatory Hilbert transform was already established 
by Wong [35]. Instead of (jl.4p . he considered H~^{f{t)){x). However, it is assumed that 
/ is a locally integrable function on [0, oo) and has an asymptotic expansion of the form 

oo 

/(t) ~ e*^*^a,i-"-", as t ^ oo, 

where < q < 1 and c is a real number. Thus, one may have c = u. Let '00(0 = f{i) 
and define V'n(i) by 

n-l 

/(t) = ^a,e*^*t-^-" + V^„(i), n>l. 

s=0 



3 



It was then shown that (see [65\ Thm. 1]) 



n-l 



H^{fm-) = EaA-)T.^-T.^ + i-M^)^ (1-6) 

for < a < 1, where 

Ea,,cix) = ^ [e-^°^r(l - a)r{a,icx) + in] , 

/•CO 

bs = / t'-^Mt)dt, s > 1, 

^0 

and 

dn[x) = + — dt, n = 0, 1,.... 

Jo 

Here T{a) is the Gamma function (21 p. 255] and r(a, z) is the complementary incomplete 
Gamma function [21 p. 260]. The expansions when a = 1 are also derived in a similar 
manner; see [65\ Thm. 2 and Thm. 3]. Moreover, the following bounds for 6n{x) are 
achieved: 

\Sn{x)\ < 

for all x > e, where M„ is a positive constant. The difficulty in applying expansions (|1.6p 
is that, as also pointed out by Wong, the coefficients bg are inconvenient for calculations. 
In a later paper [32], Ursell generalized the results of Wong and showed that these 
coefficients can be readily determined whenever the Mellin transform of f{t) is known. 
An interesting example was given for f(t) = J^^t) with applications in water waves, 
where Jo{t) is the zeroth-order Bessel function of the first kind. 

For the numerical aspects of (|1.4|) . King et al. |18j constructed a fairly robust 
numerical procedure by using convergence accelerator techniques. Unfortunately, this 
series acceleration method may suffer from difficulties when the singularity is embedded 
in a region of extreme oscillatory behavior. 

The purpose of this paper is two-fold. On the one hand, we shall derive asymptotic 
expansions of one-sided oscillatory Hilbert transforms (|1.4p as w — >• oo. To the best of 
our knowledge, none of the studies are available in this direction. Such an expansion 
clarifies the behavior of ()1.4p for large uj and also provides a powerful mean for the 
design of effective computational methods. On the other hand, in view of the fact that 
asymptotic expansions are not suitable for numerical calculation, we present efficient 
quadrature rules to approximate such integrals. It comes out that these rules depend on 
the position of x. This can be seen from (|1.4p and (|1.6p . where the integral may tend to 
zero as X — 7- oo and blow up as x — ?• 0. 

The rest of this paper is organized as follows. We perform asymptotic analysis of 
oscillatory Hilbert transforms in Section [2l The analyticity of / is of importance in our 
derivation. In Section [3l we propose efficient and affordable approaches for numerical 
evaluation of such oscillatory transforms. These methods are designed for three regimes, 
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that is, X = 0{1) orx^l, 0<x^l and x = 0, which cover all the situations. Numer- 
ical experiments show that the convergence of the proposed methods greatly improves 
when the frequency u increases. Some ideas in this paper can also be extended to study 
oscillatory Hilbert transforms with Bessel oscillators. We give a brief description of this 
aspect in Section [H We conclude this paper with some final remarks in Section [5j 

2 Asymptotic analysis of oscillatory Hilbert transforms 

2.1 Large u expansion with x > 

We start with the derivation of asymptotic expansions of oscillatory Hilbert transforms 
()1.4p for large u with x > 0. An important ingredient in our analysis is the following 
lemma which allows us to reduce the Cauchy principal integrals (jl.4p to ordinary integrals 
under certain restrictions on /. 

Lemma 2.1. Let f be a locally integrable function on [0,oo) and continuously dijjeren- 
tiable over (0,oo). Suppose that f has an analytic continuation in the first quadrant of 
the complex plane, except possibly a branch point at the origin, and there exist constants 
M > 0, 5 < 1 and < d < u such that 

\f{z)\<M\z\^e'^^'^'^'\ (2.1) 

as \z\ — )• oo in the first quadrant. Then we have 

ii-+(/(t)e^"*)(x) = e'^^'^^dt = i7re*"'V(x) + H e^'P^^dp (2.2) 
Jo t-x Jq p + ix 

for each x > 0, whenever the integral exists. 

Proof. Let us consider a quarter of the disc centered at the origin with radius R, which 
lies in the first quadrant, and denote by its boundary, i.e., 

■.= {z\\z\= R, Re (z) > 0, Im (z) > 0}. (2.3) 

For each fixed x > 0, we can find R large enough such that a half disc U^^ := {z\ \z—x\ < 
e, Im (z) > 0} can be excluded from the quarter of the disc, where e is a small positive 
number. The obtained domain is bounded by curves orientated in a counter-clockwise 
manner as illustrated in Figure [TJ 

According to our assumptions on /, it is easily seen that the integrand e*'^*|J^ of 
()1.4p is analytic in a quarter of the disc with a small indentation at x, as described above. 
We then obtain from Cauchy's theorem that 
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X — e X X + e R 

Figure 1: A quarter disc in the first quadrant with a small indentation at x. 



where dU^,, stands for the boundary of U^^^. This, together with ()1.4p . implies 

iR 



H+{f{t)e^-'){x) = hm 



+ 



dm 



t — x 



We next evaluate the three integrals on the right hand side of ()2.5p . 
A simple change of variable 



t = Re''', 0<9<-, 



vr 



yields 



t — X 



fiRe 



Re^^ - X 



Re'^idO 



< R 







Jo 


Re'^^ - X 



< 



R 



< 



R-x 
MR^+^ 



-LdR sin ( 



dO 



\f{Re'')\de 



-{u]—d)Rsin8 



(2.5) 



(2.6) 



R-x 

where in the last step we have made use of ()2.ip . Recall the well-known inequality 
sin6' > f 61, if < 6* < f ; cf. [H p. 223]. For R large enough, we obtain 



e'-^I^dt 
t — x 



< 



f\-ii.^d)Re^g^ 
Jo 

( 1 _ e-(^-'^)«') ^ 
2{u;- d){R-x)\ ; ' 



ttMR^ 



as i? — )■ 00. 



(2.7) 



It is also easily seen that 



lim e'^'^dt= lim / e-^^^^idp 

R-^oo Jq t — X R-^oo Jq ip — X 



CO 



e-^P^^dp. (2i 
p + IX 

To evaluate the third integral over the contour dU^^, we appeal to [Ml (2-7) and (2.9) 
which gives 

^^.t fit) 



lim/ e"^'^^^dt = -me'^^'^fix). (2.9) 
^^%c/+, t-x 

Finally, substituting into ([23]), we obtain ([22]). □ 



Now, we are ready to prove 

Theorem 2.2. Let f be a function as given in Lemma \2.1\ and assume that f takes an 
asymptotic expansion of the form 

oo 

/(t)~^a,t^-" (2.10) 

j=0 

as t —7- 0"*", where < a < 1. Then the one-sided oscillatory Hilbert transforms p.4p 
can be expanded in the following fashion 

e=o yj+k=e 

as u ^ oo, for each x > 0. 

Proof. By ()2.10p . it follows that, for each fixed x > 0, 

f{ip) 1 / °° \ ( °^ {ipY\ 

p + ix ^ Jx \ ^ ^^^^ I I ^ j 



oo 

-i 



j=0 / \k=0 

a,: 



i=0 \j+k=i 



oo 



£=0 \j+fc=£ 
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as p — 7- O"*". In view of (|2.2|) . an appeal to Watson's lemma (cf. |37[ p. 20]) gives us, 



oo 



1=0 \j+k=i ^ / -^0 



oo 



oo 



e=o \j+k=e 



as a; —7- oo, which is ()2.1ip . □ 

Remark 2.3. It is worth noting that the expansion (j2.1ip is only uniformly valid for x 
bounded away from 0. To clarify the behavior when < x ^ 1, we make the following 
decomposition: 

Jo t-x Jo t^{t-x) 



with fait) = t"f{t). For the first integral on the right hand side of ()2.13p . in view of 
the expansion 



fait) - faix) 



t — X 

as t 0+, then by [351 Thm. 1, p. 199] we find that 



Y,a,x^-'\ + lYa,xJ-']t + Oit' 



fait) - faix) (A ,_,\ r(l-«)ef(^-")- , X . 

-e dt^l}_^ajx^ + 1 w^oo. 



f^it -X) J I ^l^a \ ^2 

where the sum on the right side converges when < x <C 1. For the second one, we 
obtain from |35] and (13.61) below that 



oo ^iu^t r [e-^°"r(l - a)r(a, iujx) + in] , ifO<a<l, 



f^it-x) \ e*'^^[Ei(l,iwx) + i7r], if a = 0. ^ ' ^ 

Combining these results, we can see clearly that the asymptotic behaviour of the oscilla- 
tory Hilbert transforms ()1.4p depends strongly on the behaviour of the product of oj and 
X. For example, the asymptotic behaviour of (|1.4p as x — t- 0"^, a; — ?• oo and ojx — t- 0"^ can 
be derived by taking into account the asymptotic expansions of Ei(l, iujx) and r(a, iux) 
respectively. We omit the details here and only give a leading term 



fit)_^^.t^^^l ^ sd^e— - + ivr +Oiu;-), ifO<a<l 



[air) 



t-x I -aolog(a;x) + 0(1), if a = 0. 



(2.15) 
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Remark 2.4. In Theorem \2.2\ we require that / has an analytic continuation in the 
first quadrant and satisfies the growth condition ()2.ip . These conditions can be further 
relaxed to allow / has a simple pole in the first quadrant or / only has an analytic 
continuation around the real axis and satisfies (|2.ip . In both cases, the only contribution 
will be a exponentially small term in w, thus, the expansion (j2.1ip still holds. 



2.2 Large expansion with X = 0: asymptotics of Hadamard finite-part 
integrals 

When X = 0, the integrand of the one-sided oscillatory Hilbert transforms ()1.4p has a 
singularity at the origin. If the integral is divergent, the transform should be understood 
as a finite-part integral in the Hadamard sense. Note that the definition of Hadamard 
finite-part integral does not change the value of a convergent integral. It is the aim of 
this section to find the asymptotics of ()1.4p with x = 0. 

Assume that / still admits an asymptotic expansion near the origin as given in (j2.10p , 
we then formally have 



/ 

Jo 



e"^'^dt = ao —rTdt+ ^'^^i^ ^ dt. (2.16) 



There are two integrals on the right hand side of ()2.16p . The first one is divergent and 
should be interpreted as a Hadamard finite-part integral over the positive real axis. The 
integrand of the second one has an integrable singularity at the origin. Hence, it is well 
defined. 

To this end, one needs to extend standard Hadamard finite-part integrals for the 
finite interval (cf. [191 Sec. 1.4]) to semifinite integrals. Here, we adapt the definition 
from 1251: 



Definition 2.5. Let g{x) be of class C"^"'"^[0, oo) and such that 
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< CO, /c = 0, . . . ,m + 1, (2.17) 



for all p > 1. Then for any ?] > 1, a finite-part integral of order b for the positive real 
axis is defined as 

fi^W,*, (2.18) 



where |f stands for a Hadamard finite-part integral and h is an arbitrary positive constant. 

In [25], a function g satisfying ()2.17p is called allowable. The Definition 12.51 defines a 
Hadamard finite-part integral for the infinite interval from a finite-part integral on the 
finite interval, and it is independent of the choice of h. 

With Definition 12.51 we observe that 
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Lemma 2.6. 

roo „iujt 



Jo 1 -7 -logo; + if, ifa = Q, 



where 7 is the Euler constant. 

Proof. From (j2.18p . it is readily seen that 



^,dt:=f —T-rdt+ / —rdt, (2.20) 
Jo Jb 

where h is an arbitrary positive constant. If < q < 1, an appeal to integration by parts 
gives us 

b Aujt Au)e „iu]b a, , rb 

' dt = — -^ + - t-^e'^'dt 

ioj /r(i -q) 







+ 


iuj 




ab°' 


a 




^iuib 


+ 


ioj 




ab°' 


a 



ij^iJ)^ - -i^h) - t-"e*-*dt) , (2.21) 



and 



= — + — / t-^e'^^'dt = — + — 6^-"Ei(a, -zw6), (2.22) 



ab"" a Ji, ab'^ a 

where e is a small positive number and 

/oo 
t~Pe-'^dt, p>0, Re(z)>0, (2.23) 

is the exponential integral. In (|2.2ip . by neglecting the divergent term and noting 
that the last integral vanishes as e — t- 0"*^, we obtain 



Combining ([2:20]) . and we get the desired result. 

Similarly, if a = 0, we note that 

f^b rooX 'lujt 



j +j J — dt = Ei(l,-ia;e) = -7-loga; + i--loge + 0(e), (2.25) 



as e — )• 0, where 7 is the Euler constant. Hence 

foo Aujt 



^ := -7-loga; + i-, (2.26) 

JO ^ 

as shown in (j2.19p . □ 
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A combination of Lemma 12.61 and (|2.16|) gives us 

r = / ^"'T"" r(l - a)ao + f,^ e-*iW^dt, if < a < 1, 

Jo t \ (-7 - logL^ + if )/(0) + e'^^' ^^'^'/^^U t, ifa = 0. 

(2.27) 

Here we have made use of the fact that ao = /(O) if a = in (|2.10|) . To derive 
asymptotic expansions of Hadamard finite part integrals (j2.27|) . it is then sufficient to 
find asymptotics of Fourier-type integrals f(t)e^^^dt, where 

:= fit)-^(^ot~'" ^ < a < 1. (2.28) 

Our result is stated below: 

Theorem 2.7. Let f be a locally integrable function on [0, oo) andm times continuously 
differentiahle over (0, oo), m being a positive integer. Suppose that f has the asymptotic 
expansion ()2.10p . and this expansion can be differentiated m times. Moreover, each of 
the integrals 

/oo 
f^'\t)e'^^dt, s = 0,l,---,m, (2.29) 

converges uniformly for u sufficiently large, where f is given in (j2.28p . Then, as u ^ oo, 
we have 



i 



^i^^r(l - a)ao + E a,ef + o(a;-), z/0 < a < 1, 

(if -7-logw)/(0) + E«^ef^^ + o(w-'"), ifa = 0. 



(2.30) 



Proof. From ()2.28p and ()2.10p . it is readily seen that 



/» = M^^~j:a,+,t^-" (2.31) 

as t — )• O"*". This, together with [351 Thm. 1, p. 199], leads us to the following asymptotics 
of Fourier integrals: 

r me^^' = E + o(^-). (2.32) 

The asymptotic expansions (j2.30p of Hadamard finite-part integrals then follows from a 
combination of (p?^ and ([232]) . □ 

Remark 2.8. We do not require that / has an analytic continuation to the first quadrant 
in Theorem 12.71 However, such requirement is essential in the proof of Theorem 12.21 
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2.3 Some examples 

We conclude this section with appHcations of Theorems 12.21 and 12.71 to some concrete 
examples. 

Example 2.9. Consider the integral 

-^°°e*^*|-^dt, c>0. (2.33) 

This is the one-sided oscillatory Hilbert transform of f{t) = e~'^*. Clearly, / is an entire 
function and satisfies all the assumptions of our theorems. Since 

fit) - E (2.34) 



as 



t — )• 0"*", one has 



a = 0, aj=^-4^, j = 0,l,2,---, (2.35) 



in the notation of ()2.10p . Hence, from ()2.1ip and (|2.30p . it is easily seen that 



Jo t-^ t". L.4i 



£=0 \j+k=£ •' 

and 



x>0, (2.36) 



f e-^'-^dt ~ - 7 - log - + E ^e^^^A' (2-37) 

Jo t 2 t:^ ^ ^ 



as — 7- oo. In particular, if c = 0, i.e., / = 1, the above two expansions can be simplified, 
and we have 

Jo t-x {u:xY+^ ^ ' 



i=0 

and 



roc giwt ^ 

J. dt = i---f-loguj. (2.39) 

Jo ^ 

It is worthwhile to point out that the asymptotic expansion ()2.38p is the same as the 
large x expansion of jr^dt with w > 0; see [32l Lem. 1]. 

Example 2.10. As the second example, we consider 

/•oo /7 

t e''^*7^ TT ^dt. (2.40) 

Jo il + t){t-x) ^ ^ 
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Thus, fit) = Vi/{1 + t), which can be analytically extended to the complex plane with 
a pole at —1 and a branch cut along the negative axis. As t — )• 0"'" , we have 

oo oo 

fit) ~ VtY^i-t)^ = (2.41) 

j=0 j=l 

Hence, in the notation of (|2.10p . a = 1/2, ao = 0, aj = (— 1)-'"^^, j = 1, 2, • • • . An appeal 
to ([2TT1) and (fOO]) then gives 

Jo ii + t)it-x) i + x ^ jH yf^^ xi j 

for X > 0, and 

Jo til+t) ' J-h 

3 Computation of oscillatory Hilbert transforms 

Although the asymptotic expansions derived in the previous section provide essential 
insights into the behaviour of the oscillatory Hilbert transforms for large lo, they are not 
suitable for computational purpose, since asymptotic expansions are typically divergent, 
one can not simply keep on adding terms of the expansion in order to improve the accu- 
racy of the approximation. In this section we shall focus on fast numerical computation 
of oscillatory Hilbert transforms (|1.4p . According to the position of x, we classify our 
discussion into three regimes, namely, x = 0(l)or2;^l, 0<a;^l and x = 0. 
Since these regimes exhibit different asymptotic behaviour, they also require different 
numerical methods. 

3.1 The regime x = or x ^ 1 

If X is not so close to the origin, say, x = 0(1) or x ^ 1, the one-sided oscillatory Hilbert 
transform p.4|) can be approximated efficiently using the generalized Gauss-Laguerre 
quadrature rule. To see this, we observe from ()2.2p that 

H+ifit)e'^^)ix) = mfix)e''^'' + - / e"" ! dq 

^ Jo Z^+^x 

= -/(x)e-^ + '-^^^t3pl r q~-e-^^dq, (3.1) 

with fait) = f^fit) and where a is defined as in ()2.10p . We have made use of the 
change of variable q = up in the first equality. On account of (j2.10p . it is easily seen 
that fai^)/i^ + ix) behaves like a polynomial near the origin, if x = 0(1) or x ^ 1. 
This is exactly the situation that can be handled by the generalized Gauss-Laguerre 
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quadrature rule. Let {tj,Wj}'j^^ be the nodes and weights of the generahzed Gauss- 
Laguerre quadrature rule with respect to the weight t~"e~*, with < a < 1. Then the 
one-sided oscillatory Hilbert transform H'^{f{t)e^^^){x) is approximated by 

Qn{f,uj,x) = nTf{x)e + -— — > Wfc , , . (3.2) 

^ — + IX 

k=l 

Applying the error expression of the n-point generalized Gauss-Laguerre quadrature rule 
|10^ p. 223], we can estimate, for each fixed x, the quadrature error as follows: 




exp(-f TTf) n!r(n - a + 1) ( fo.{i^ 



(2n)! \3- + ix 



' ^ + ix ' ^ +ix 

" ^ k=l <^ / 

(2n) 



9=€ 



= 0(a;-^"-^+") (3.3) 

for some constant € C as w — ?• oo. Note that one gains the factor u]~'^^ from taking 
the (2?i)-th order derivative of a function of q that depends only on q/oo. We observe 
that the accuracy of the quadrature rule (3.2) rapidly improves as uj grows. This result 
may come as a surprise here, since we started out with a highly oscillatory integral that 
typically requires a growing number of quadrature points to evaluate as the frequency 
increases. Yet, we like to point out that this result is entirely parallel to the case of 
finite Fourier integrals of the form (jl.ip , where an application of Gauss-Laguerre yields 
similar behaviour for large oo. In what follows, we give several examples to illustrate 
the convergence of our quadrature rule Qnif,u},x). Throughout this paper, all the 
computations have been performed using Maple 14 with 32-digit arithmetiqj. 

Example 3.1. Let us consider ()2.33p with c = 1, that is, f{t) = e~*, and one has 
a = in (j3.2p . The exact solution of ()2.33p can be obtained by using the definition of 
the Cauchy principal value integral and taking the series expansion of the exponential 
function. The result is 



i 



oo p—ct 
^iu^t^ 

t — X 



Ei(l, (c - ^.)x) - 2 X; E ^' ( , J ^ + t) 

^s(-) - ^ E (E K'^r ) 



(3.4) 



where 



Si(x)= r^dt, x>0, (3.5) 
Jo * 



*The use of increased precision is just to show the convergence rates of our methods. 
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is the sine integral; cf. [5]. In the case of c = 0, one can check that 



i 



t — X 



dt = e*""^ (Ei(l, -iuox) + i2Si(a;x)) = e*^^' (ivr + Ei(l, iwx)) . (3.6) 



We compute the error for x = 1 and x = 5 with different frequency oj and n ranging 
from 2 to 10. The results are illustrated in Figure [2j 



2 3 4 5 6 7 



° co=16 • 0)^8 * (0=4 ° co=2 



2345 6 789 10 

° co=16 • a>=8 * co=4 ° a=2 



Figure 2: The error of the quadrature a;,x) for x = 1 (left) and x = 5 (right) with 

f{t) = e~* and n ranging from 2 to 10. 



Example 3.2. We next consider the function f{t) = which grows exponentially in 



the complex plane. We can check that it satisfies ()2.ip with d = 1, and 0=3 
The error is illustrated in Figure EJ 



in M. 



(B=16 • a)=8 * co=4 ° (B=2 



° (0=16 • (0=8 * (0=4 ° (0=2 



Figure 3: The error of the quadrature (5„(/, a;,x) for x = 1 (left) and x = 5 (right) with 
f(t) = and n ranging from 2 to 10. 
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° 05=16 • (0^8 * £0=4 ° co=2 



° a)=16 • 05=8 * 05=4 ° (xs=2 



Figure 4: The error of the quadrature Qn{f,u},x) for x = 1 (left) and x = 5 (right) with 
/(O = j+t ^ ranging from 2 to 10. 

Example 3.3. We return to Example 12.101 i-S-, f{t) = V^/(l + t), which corresponds 
to a = 1/2. Numerical results are displayed in Figured! 

All these examples show that the accuracy of quadrature rule (|3.2p improves rapidly 
as n increases. Meanwhile, the convergence is faster for larger uj. 

3.2 The regime < x < 1 

When X is close the origin, i.e., < x ^ 1, the accuracy of the generalized Gauss- 
Laguerre rule deteriorates since the integrand on the right hand side of ()3.ip is nearly 
singular. To this end, we make the following decomposition of oscillatory Hilbert trans- 
forms: 

H+{f{t)e''^'){x) = r ^^e''''dt+ r I^e'^'dt, (3.7) 
Jo 't — X t — x 

where a is a positive number larger than x. Let the two integrals on the right hand side 
of (|3.7p be denoted by Ii{x) and hix), respectively. Note that the integral h^x) is no 
longer singular. A fast computation of oscillatory Hilbert transforms in this regime is 
then reduced to the numerical study of Ii{x) and l2{x), which will be discussed in the 
next two sections. 
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3.2.1 Computation of /i(x) 

The integral Ii{x) is a finite oscillatory Hilbert transform. By scaling the interval of 
integration to ( — 1, 1), we have 



(x) = f 
Jo 



t — X 



-1 y-T 

1 {y+l)-^h{y) 



1 



y-T 



where a; = ^f,r = ^- lG (-1, 1), and 

h{y) = fai^iy + 1)) = + 1))" / {^{y + i] 



(3.8) 



(3.9) 



The computation of ()3.8p with a = has been discussed in [33]. Let PNiy) be 
the Lagrange polynomial which interpolates h{y) at the Clenshaw-Curtis points yj = 



cos 



N ;> 



j = 0, . . . , A^. Then one has (see [6]) 

TV 



(3.10) 



fc=o 



with 



AT 



N 



Y,"Hyj)Tj{yk), k = o,...,N, 



(3.11) 



j=0 



where the double prime denotes a sum whose first and last terms are halved and Tj{y) 
is the Chebyshev polynomial of the first kind of degree j. The coefficients can be 
computed efficiently by FFT [11]. Replacing h{y) in p.Sp by PAr(y), we have 



' fjy) 
ly-T 



e'^^dy 



^ PN{y) 



e'^ydy 



■ 1 y-T 
^ PN{y) -pn{t] 



y-T 
rk-i 



e'^^dy+pNir) 



1 ^iCuy 



fe=o 



.n=0 



-ly-T 

•1 

+ Pn{t 



dy 

^dy, 
^ly-r 



(3.12) 



where 



Mn = j^Un{y)t 



e^^^dy, n > 0, 
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and where Un{y) is the Chebyshev polynomial of the second kind of degree n. The last 
integral in ()3.12p can be computed in closed form and the M„ can be computed by using 
a three-term recurrence relation. The advantage of this method is that it converges 
uniformly with respect to r G (—1, 1) as — t- oo, provided h{y) is analytic in a small 
neighborhood containing [—1,1]. However, this method costs 0{N'^) operations which 
make it inefficient for large N . 

In the following we shall present a more efficient implementation, which costs only 
0{N \og2^) operations, to compute the finite oscillatory Hilbert transform (jS.Sp . The 
key observation is that we can write {pN{y) ~ Pn{t)) / {y — t) in terms of Tk{y) and this 
process can be performed in only 0{N) operations. 

We approximate Ii{x) by 

g«(/,^,x) = (^)'%- (£(y + i)-"^e-^d,) 

-Y^e^H tiy + iy-^-^^^^^^^^e^'^ydy 
2/ \J^i y-T 

(y + 1)-- dy). (3.13) 

I y-r J 

We next expand {pN{y) — Pn{t)) / {y — t) in terms of Tfc(y) and obtain 

'""^'^-''^^^ = Y.'b^Myl (3.14) 
y — T ^-^ 

where the prime denotes the summation whose first term is halved. The coefficients 
satisfy a three-term recurrence relation 

6f_i = 2af + 2r6f A; = iV - 1, . . . , 1, (3.15) 

and the first two initial values are given hy = and b^_i — a^; see |13j . Inserting 
(I3TI]1 into (f3lB gives us 

Q«(/,L^,x) = (-) e- j;'&^zf)+p^(T)-/ (y + 1)- dy], (3.16) 

where 

4°) := Zk{a,Cj) = j^{y + l)-"Tfe(y)e^^J'dy, A; > 0. (3.17) 

Now, we only need to compute the moments and the integral on right hand side 

of ()3.16p efficiently. It turns out that Z^""* satisfy the following four-term recurrence 
relation (see [31]): 

iOj{n - l)zi"\ + [2(n - Q + l)(n - 1) + iuj{n - 2)]^") 

+ [2n(n + a - 2) - iu:{n + l)]4-i - i^nZ''n-2 = -2^-'^e'^ (3.18) 
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with the first three initial values given by 



Zq - 




~i 








Z2 - 







where 



e5-^(i-°)7(l-a,-2ic:.), 
el-(2-°)^(2 - a, -2ia}) - 
el-(3-")^(3 - a, -2iCo) - ^zf - Szj"^ 



7(a,z)=/ t"~^e"*dt, Re(a)>0, (3.19) 



is the incomplete gamma function [2], p. 260]. The recurrence relation ()3.18p can be used 
to calculate Z^"^ stably in the forward direction provided N < 2Cj. If > 2(1;, we can 

compute the additional moments Z^"'' , 2uo < k < N stably by solving a boundary value 
problem with two starting values and one ending value [31]. To see this, let uq := [2a)], 
where [ ] denotes the integer part, and choose a positive integer A^i > maxjno, A'^}. We 
then define a matrix A = {aj^k)flZi'^~^^ of size A'^i — no + 1 by 

iu;{no+j - 2), k = j + l, 

_i 2(no - a + j)(no + i - 2) +zLD(no +i - 3), k = j, , . 

Hk \ 2{no+j-l){no + a + j-3)-iu{no+j), k = j - 1, ^''•^"^ 

-ioj{no+ j - 1), k = j-2, 

and a vector b = (61, 62, • • • , &Afi-no+i)^ by 

_22-"e^-_(2no(no + a-2)-ici(no + l))zi"Li + icino4"L2, i = 1, 



.22-"e- + iLD(no + l)z£Li, j = 2, 



-2^-"e*-, i = 3,...,iVi 



22-"e^- - iu{N, - l)4^Vi' j = Ni- 710 + 1 



(3.21) 

where the superscript stands for transpose. Here, the moments Z^^_2 and Z^^_^ can 

be obtained by using the forward recurrence (j3.18p . while the value of -^^^+1 hi the last 
component can be set equal to zero if the selected parameter A'^i is sufficiently large. We 
have that the vector consisting of the additional moments defined by 

^ - l^no ' ^no+1' • • • ' ^ATi ) ' 

is the unique solution of the linear system 

= b. (3.22) 
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To compute the integral on right hand side of (|3.16|) . we note that 



(y + ly^ dy = 7T— ^(^vr + e-*-'^r(l - a)r(a, icD(l + r))) 

1 y-T (l + r)" 

I 

, e"* ^ r^dt, (3.23) 

^ Jo (2 + |)-(l-T + f) ^ ^ 



where the last integral is free from singularity as r — )• —1, and can therefore be con- 
veniently evaluated by Gauss-Laguerre quadrature rule. For the case a = 0, it can be 
further written in a closed form (see [5]): 



1 giiliy 



dy = cos(a;r) [Ci(tii) — Ci('U2)] — sin(LDr)[Si(ui) + Si(n2)] 



ly-r 

+ i [sin(wT)[Ci(ni) - Ci(n2)] + cos(wr)[Si(ni) + Si{u2)]] , (3.24) 



where 

Ci(x) = - / ^dt, x>0, (3.25) 

Jx t 

is the cosine integral, ui = — t) and tt2 = a;(l + r). 

Based on the above discussion, we outline our algorithm with estimation of compu- 
tational complexity as follows: 

Algorithm I: Computation of h{x) 

1. Choose a positive constant a larger than x and A'^i > maxjrio, A^}. 

2. Compute the coefficients {a^}k=o in (l3TT]l by using FFT with C'(A^log2 A^) oper- 



-'fc ik=0 

ations. 



3. If < 2lo, evaluate the moments {zj^^}^^^ from the recurrence relation ()3.18p 

''k Jk=no 



in 0{N) operations. If > 2Cb, we compute the additional moments {Z^"^}^ ^ 



by solving ()3.22p . in 0{Ni) operation^. 

4. Compute Pn{t) from its barycentric form in 0{N) operations [1]. 

5. Calculate Q%\f, Cj^x) in 0{N) operations by the Clenshaw algorithm: 

5_i = 0,5o = ^4"^^o = 0, 



Wk = Wk-i + 2a^Sk^^, fc = l,2,...,Af-l, 

'1 iQy 



— ^fc"^ + 2TSk-l — Sk-2, 



^This can be done with a minor adaptation of Oliver's method for the LU decomposition of a tridi- 
agonal matrix [23 ■ We omit the details. 
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The integral + 1) "'^z^dy can be evaluated efficiently from (|3.23p by using 

a Gauss-Laguerre quadrature rule. 

The total computational complexity of computing /i(x) is C'(A^log2iV) operations 
if < 2a) and is C'(A^log2 A^) + 0{Ni) operations if > 2(1;. Therefore, even when 
we choose A'^i = C'(A^log2 N), the total computational complexity of computing Ii{x) is 
still 0(A^log2 N). Meanwhile, this cost is independent of u. 

In the case a = 0, i.e., / is analytic in the neighborhood of the origin, the above 
algorithm can be further simplified. Indeed, if a = 0, we have 

= 4 fe^" - (-l)^e-'^) - ^Mfc.i, k>l. (3.26) 

Since Mj. satisfies a three-term recurrence relation (see [HJ [33] ) 

21 2 / \ 

Mz + T^Mi„i-M,_2 = T^(e^^-(-iye^'^), l>2, (3.27) 

lUJ lUJ \ / 

with initial values Mn = Mi = 4i(^ - the four-term recurrence relation 

(j3.18p can readily be reduced to a three-term recurrence relation. In practice, one can 
evaluate from (|3.27p for N < oj. If N > Cj, the additional moments M^, u < k < N 
can be evaluated stably by solving a tridiagonal system [8] . We then have the following 
simpler algorithm: 

Algorithm II: Computation of Ii{x) with q = 

1. Choose a positive constant a larger than x. 

2. Compute the coefficients {a^^l^o using FFT in C'(A^log2 A^) operations. 

3. Evaluate the moments {M}S\^~^ by (|3.27p li N < uj. li N > uj, we compute the 
additional moments {M^j^lj by the second phase algorithm in [8]. Then we 

evaluate {4°^}fcL~o^ dSSi), in 0{N) operations. 

4. Compute Pn{t) from its barycentric form in 0{N) operations 1^. 

5. Calculate Q^]^\f,Cd,x) in 0{N) operations by the Clenshaw algorithm: 

5-i = 0,5o = ^4°^W^o = 0, 

Wk = Wk-i + 2a^Sk-i, = 1, 2, . . . , AT - 1, 

Q'^^\f,Cj,x) = e'^ (Wn^i + SN^ia^+PN{T)-j- ^ J—^dv) ■ 

The integral y^'^y evaluated from (j3.24p . 

The computational complexity is C'(A^log2 A^) operations. 
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3.2.2 Computation of hix) 



The computation of l2{x) is relatively easy, since the integrand does not have singularity. 
Under the same assumptions as in Lemma 1 2. 11 it is readily seen that 



hix) 



t — X 



f{a + ip) ^ 
a — X + ip 



dp. 



(3.28) 



As a is an arbitrary real number larger than x, we may select it so that the integrand of 
the last integral in ()3.28p is well behaved. Thus, this integral can be computed efficiently 
by Gauss-Laguerre quadrature rule: 



l2ix)c^Q^^\f,oj,x) 



k=l 



(3.29) 



where {tk,Wk} are the nodes and weights of the Gauss-Laguerre quadrature associated 
with the weight e~*. 

(2) 

Example 3.4. We apply quadrature rule Qn {f,i^,x) in (j3.29p to calculate l2{x) with 
f(t) = Vi/{l + t) and x = 0.02. The error is presented in Figure [5l We can see that more 
accurate approximations are obtained as a increases. For each fixed a, more accurate 
approximations are obtained as to increases. 



5 6 7 



(0=20 • OFlO * CIfS ° OFl 



(0=20 • KFlO * KfS ° 03=1 



(0=20 • OFlO * cifS 



Figure 5: The error of the quadrature Qh, {f,u;,x) for a = 1 (left), a = 2 (middle) 
and a = 4 (right) with /(t) = Vt/(1 + t) and n ranging from 2 to 10. Here we choose 
X = 0.02. 



3.2.3 Numerical examples 

We present in this section some numerical experiments to illustrate the efficiency of 
numerical methods presented in the above two sections. 

Example 3.5. Let us consider the integral (j2.33p with c = 1 and x close to zero. The 
exact solution is given in ()3.4p . Since a = in this case, we use Algorithm II to compute 
Ii{x), and the Gauss-Laguerre quadrature rule Qn (/, w, x) defined in (j3.29p to compute 
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I2{x). The absolute error for several values of x is presented in Table 1, which indicates 
the proposed method is uniformly accurate as x — )• 0. The absolute error for several 
values of oj with a = 1 and x = 0.02 is presented in Table 2. We can see that the 
convergence is quite rapid as N and n increase. 

Table 1: Absolute error in computing (|2.33p with c = 1, a = 1, = 10 and x = 10^^. 



n 


N 


6 = 1 




6 = 2 




5 = 3 




5 = A 




4 


4 


1.22x10" 


-5 


3.80x10" 


■5 


3.52x10" 


5 


3.42x10" 


5 




8 


3.30x10-^ 


-7 


1.83x10" 


-7 


1.73x10" 


-7 


1.72x10" 


■7 




16 


3.30x10^^ 


-7 


1.83x10" 


■7 


1.73x10" 


-7 


1.72x10" 


■7 


8 


4 


1.19x10" 


-5 


3.80x10" 


■5 


3.52x10" 


-5 


3.41x10" 


■5 




8 


2.69x10" 


10 


1.09x10" 


-10 


1.04x10" 


10 


1.03x10" 


10 




16 


2.86x10" 


10 


1.09x10" 


10 


9.96x10" 


11 


9.87x10" 


11 


16 


4 


1.19x10" 


-5 


3.80x10" 


■5 


3.52x10" 


■5 


3.41x10" 


■5 




8 


2.37x10" 


■11 


2.65x10" 


■12 


9.36x10" 


■12 


8.00x10" 


■12 




16 


1.30x10" 


■14 


2.97x10" 


■15 


2.58x10" 


■15 


2.54x10" 


■15 



Table 2: Absolute error in computing ()2.33p with c = 1, a = 1, x = 0.02 and several 
values of oj. 



n 


N 


to) = 5 




o; = 20 




= 80 




uj = 320 




4 


4 


2.84x10" 


■5 


2.30x10" 


■5 


1.56x10" 


5 


1.73x10" 


-5 




8 


1.81x10" 


■5 


8.92x10" 


10 


1.28x10" 


11 


1.81x10" 


-11 




16 


1.81x10" 


5 


8.69x10" 


10 


4.89x10" 


15 


1.92x10" 


-20 


8 


4 


1.65x10" 


5 


2.30x10" 


5 


1.56x10" 


■5 


1.73x10" 


-5 




8 


8.16x10" 


■10 


3.12x10" 


■11 


1.28x10" 


11 


1.81x10" 


-11 




16 


1.08x10" 


10 


2.00x10" 


■14 


8.08x10" 


■24 


3.65x10" 


-25 


16 


4 


1.66x10" 


5 


2.30x10" 


5 


1.56x10" 


-5 


1.73x10" 


-5 




8 


8.69x10" 


11 


3.12x10" 


11 


1.28x10" 


11 


1.81x10" 


-11 




16 


7.49x10" 


11 


5.44x10" 


■21 


2.78x10" 


■25 


3.65x10" 


-25 



Example 3.6. We consider (j2.40p with x close to the origin. Since a = 1/2 in this 
case, we use Algorithm I to compute Ii{x) and the Gauss-Laguerre quadrature rule 
Qn (f,uj,x) in ()3.29p to compute l2{x). For simplicity, we choose A^^i = 2A^ when 
> 2a; in our implementation and the last integral in (|3.23p is computed by using 32- 
point Gauss-Laguerre quadrature. The absolute error is presented in Table 3 for a = 1 
and u = 10, which implies that the convergence is quite rapid as N and n increase. In 
Table 4 we present the absolute error for several values of oj with fixed a and x. As we 
can see, the accuracy greatly improves as uj increases. 
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Table 3: Absolute error in computing ()2.40p with a = 1, u = 10 and x = 10 



n 


N 


5 = 1 




5 = 2 




5 = 3 




6 = 4 




4 


4 


1.53x10-^ 


-3 


3.56x10" 


-3 


4.56x10" 


-3 


4.67x10" 


-3 




8 


1.31x10-^ 


-6 


1.63x10" 


-7 


2.87x10" 


-6 


3.25x10" 


-6 




16 


1.16x10" 


-7 


4.51x10" 


-8 


4.08x10" 


-8 


4.03x10" 


-8 


8 


4 


1.53x10" 


-3 


3.56x10" 


-3 


4.56x10" 


-3 


4.67x10" 


-3 




8 


1.23x10" 


-6 


1.29x10" 


-7 


2.89x10" 


-6 


3.27x10" 


-6 




16 


9.20x10" 


-11 


2.22x10" 


-11 


1.62x10" 


-11 


1.51x10" 


-11 


16 


4 


1.53x10" 


-3 


3.56x10" 


-3 


4.56x10" 


-3 


4.67x10" 


-3 




8 


1.22x10" 


-6 


1.29x10" 


-7 


2.89x10" 


-6 


3.27x10" 


-6 




16 


1.39x10" 


-12 


2.46x10" 


-12 


1.33x10" 


-12 


2.38x10" 


-12 



Table 4: Absolute error in computing ()2.40|) with a = 1, x = 0.02 and several values of 



n 


N 


uj = 5 




o; = 20 




W = 80 




UJ = 320 




4 


8 


5.50x10"^ 


6 


2.02x10" 


6 


2.04x10" 


6 


2.38x10" 


■6 




16 


5.08x10"- 


6 


2.15x10" 


-10 


1.23x10" 


12 


2.60x10" 


■12 




32 


5.08x10- 


6 


2.15x10" 


■10 


1.19x10" 


■15 


4.65x10" 


■21 


8 


8 


2.06x10" 


'6 


2.02x10" 


-6 


2.04x10" 


■6 


2.38x10" 


6 




16 


2.41x10" 


-8 


4.17x10" 


■13 


1.23x10" 


■12 


2.60x10" 


■12 




32 


2.41x10" 


■8 


3.86x10" 


■15 


2.59x10" 


■24 


6.75x10" 


■25 


16 


8 


2.08x10" 


6 


2.02x10" 


■6 


2.04x10" 


■6 


2.38x10" 


■6 




16 


1.37x10" 


11 


4.13x10" 


■13 


1.23x10" 


■12 


2.60x10" 


■12 




32 


1.39x10" 


■11 


9.01x10" 


■22 


1.11x10" 


■24 


6.75x10" 


■25 



Remark 3.7. When a = 0, the uniform convergence of the quadrature rule ()3.13p has 
been proved in [SB]. When < a < 1, however, we don't currently know whether the 
uniform convergence still holds. Numerical results show that the quadrature rule ()3.13p 
remains accurate when x is small, as can be observed from Table 3. 



3.3 The regime x = 

When X = 0, we need to deal with the Hadamard finite-part integral e^'^^^Y^dt, 
which is introduced in Section [2.21 In view of ()2.27p . it suffices to find a fast method for 
evaluating the integral 



f{t) - aot 



-a 



t 



e"^'dt, (3.30) 
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or equivalently, 

/•oo 

t-'^g{t)e'^^dt, (3.31) 



where 

Since g{t) is holomorphic in the first quadrant, which follows from the assumptions in 
Lemma 12.11 and Theorem 12. 2^ we can apply Gaussian quadrature rules as proposed by 
Wong [36] : 

j\-g{t)e^-'dt = ^"P((^ +E^{a), (3.33) 

where {tk-,Wj}^=i are the nodes and weights of the generalized Gauss-Laguerre quadra- 
ture with respect to weight t~"e~*. The reminder En{g) is given by 

Therefore, we can approximate e^^^^^dt by 

( ^^2-a)^^a exp((l-a)f)" fitk 



ni-a)a, + '-^^^^^Y.^,g(^], .f < « < 1, 

i— 1 ^ ^ 

(3.34) 



(i--7-logt^)/(0) + ^u;fc^^^ , ifa = 0, 



u 1 
fe=l 



Example 3.8. We use quadrature rule (j3.34p to approximate (j2.27p with f{t) = e~* 
and f{t) = Vi/{1 + t), which corresponds to a = and a = ^ respectively, for several 
values of u. The absolute errors are shown in Figure [6l Clearly, the convergence of the 
quadrature rule p.34|) is rapid and satisfactory. 



4 Extensions 



Oscillatory Hilbert transforms with Bessel type oscillators are defined by 

H+{f{t)Mu;t)){x) := P^J,{ujt)dt 
Jo t-x 

and 

H+{f{t)Y,{ut)){x) := r p^Y,{u;t)dt, 
Jo t — X 

where J^it) and Y,^{t) are the Bessel functions of the first and second kind, respectively. 
Such kind of transformations have applications in physics such as water-wave radiation 
problem [23] . It turns out that Lemma 12.11 can also be extended to oscillatory Bessel 
Hilbert transforms, by using similar ideas. 
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n ts>=20 • £0=10 * CL)=5 o (0=1 



Figure 6: The error of the quadrature (fOil) for f{t) = e"* (left) and f{t) = + 1) 

(right) with n ranging from 2 to 10. 



Lemma 4.1. Suppose that f has an analytic continuation to the right half plane, except 
possibly a branch point at the origin, and there exist constants M > 0, (5 < | and 
< d < oj such that 

\f{z)\<M\z\^e'^^'^^'\ (4.1) 



as \z\ CO in the right-half plane. Then, 



and 



H+{f{t)J,{uot)){x) = -7Tf{x)Y,iux) 



H+{f{t)Y,{ojt)){x) = 7Tf{x).Mux) + 



Ku{ujy) 



I 2 , 2 92{y)dy, 

vr Jq y + 



(4.2) 



(4.3) 



whenever the integrals exist, where Ky{t) is the modified Bessel function of the second 
kind and 

9j{y) = (^ + iy)e^^"7(%) + i-iy+H^ - iy)e'^'''f{-iy), j = l, 2. 

Proof. We only give a sketched proof of ()4.2p . since the proof of ()4.3p can be handled in 
a similar manner. On account of the identity (see [21 Eq. 9.6.4]) 



iTrMz) = e-'i^'K^ize-'^^') - e'^^' K^{ze^^'), 



I are z\ < 



TT 



2' 



it follows 



f 

Jo 



fit) 



t — X 



Ju{oJt)dt 



1 

ivr 



e-f-f^ Ii^K,^_,^t)dt-e^-^-[ 



fit) 



Jo t — X 



Jo t-x 



Ki,{iujt)dt 



(4.4) 
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For the first integral on the right hand side of (|4.4|) . by considering the same contour as 
shown in Figured! we obtain from Cauchy's theorem that 



Jo 



m 

t — X 



Ky[—iu}t)dt 



ii: f {x)Ki,{—iujx) 



Ky{—iojz)dz + i 



R 



fjiy) 
iy -X 



K,{ojy)dy. (4.5) 



Recall the asymptotic expansion of the modified Bessel function K^{z) O P- 378] 



1 + 



Au' - 1 



+ 0{z- 



z — )• oo, I arg z\ < 



37T 



We have the following estimation of the integral over for large R 



m 

Y Z — X 



■Ku{—ibjz)dz 



< 
< 



Re^^ - X 



R-x 
MR^+^ HF 



R-x \j 2ioR 
0, as i? — )■ oo. 



K^i-iujRe'yRe'Ue 

^~(u^-d)Rsme^-^ + 0{R^^))de 



Hence, letting i? — )■ oo in (j4.5p . we arrive at 
f{t) 



f 

Jo 



t — X 



Ky{—iijjt)dt = it: f {x)Ky{—iu}x) + i 



I^K,{ujy)dy. (4.6) 

^y-x 



Similarly, we can deform the integration path to the negative imaginary axis for the 
second integral on the right hand side of ()4.4p and get 



Jo 



fit) 



Ky{iujt)dt = —in f {x)Ky{iujx) + i 



10 t — x 

A combination of (ji^ . (ji^ and (jiTT]) gives us (fOj) . 



fi-iy) 

iy + x 



K,{ojy)dy. (4.7) 



□ 



From the above lemma, we establish several interesting identities with z/ = 0, 1, which 
have not been found in classical reference books [21 [T^] and which might have important 
applications in practice: 



Corollary 4.2. We have 

lO 

and 



f 

Jo 



IT 



dt = i-iy+'-[U_,iu;x) + i-irY,iujx)], 1^ = 0,1, 



f 

Jo 



-dt = ■kJq(UJX) D-l^Q[pOX), 



t — X 



(4.8) 



(4.9) 
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and 

^^^-x '^^ " 7rxJi(a;x) - ^S_2,i(a;x), (4.10) 

where Yly[z) is the Struve function and S^^y{z) is the Lommel function of the second 
kind. 

Proof With / = 1 and = 0, 1 in (02]), it follows that 

M^t) _ f -nYoicvx) - f ^dy, if = 

yKi (ui 

y'^+x 



/•oo 

Jo 



t-x^'-\ -vry,(a;x) - f ^d,, if . = 1. ^'-''^ 



Recall the identity (see [la Eq. 6.566.3]) 



/ o o dy = - -. — -[YL_^{ujx)-Y_^{ujx)l cj > 0, Rex > 0, Keu > --. 

Jq + x^ 4cos(z^7r) 2 



Formula (|4.8p then follows from inserting the above identity with v = and z/ = 1, 
respectively. 

To show (|4.9p . we set / = 1 and = in ()4.3p and obtain 



JO 



— dt = ■kJq{ojx) / Trdy. (4.12) 

Iq t-x IT Jo y^ + x^ 

Since 

POO 

/ yi+'^(y2 + x^rK,{u;y)dy = 2-T{v + l)x-+^+^u'^-^ S^.,^^+,+i{ujx) (4.13) 

JO 

for Rex > 0, RetJ > and Rei^ > —1 (see [121 Eq. 6.565.7]), substituting the above 
identity with z^ = and /i = — 1 into (|4.12p gives us (j4.9p . 

Finally, the proof of ()4.10p is similar to that of ()4.9p . and we omit the details here. □ 

We expect that Lemma |4. II might play an important role in the asymptotic and nu- 
merical study of oscillatory Hilbert transforms. Note that the modified Bessel functions 
Ky have a non-integrable singularity at for u >1. 

Finally, we like to recall that the numerical method proposed in [3] generalized steep- 
est descent-based methods for Fourier-type integrals of the form (jl.ip to oscillatory 
transforms of the form 

[•oo 

f{x)H{ujx)dx, 



where H can be, e.g., a Bessel function. We expect that this method can be extended 
in turn to compute oscillatory Hilbert transforms with more general oscillators as well. 
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5 Concluding remarks 

In this paper, we have considered asymptotic expansions and fast computation of the 
oscillatory Hilbert transforms (|1.4p . Unlike previous work, which focused on the be- 
haviour of oscillatory Hilbert transforms for large x, we derive asymptotic expansions of 
such transforms for large to. These expansions clarify the asymptotic behaviour for large 
CO and provide a powerful approach for designing efficient and accurate approximation 
methods. 

Numerical methods for the calculation of the oscillatory Hilbert transforms are pre- 
sented. We classify our discussion into three regimes, namely, x = C(l) or x ^ 1, 
< X ^ 1 and x = 0. For each regime, we have designed efficient numerical approaches. 
Even for small values of w, our proposed methods remain quite accurate. Numerical 
examples are provided to confirm our analysis. 

In the implementation of our methods, the function f{x) is required to be analytic in 
the first quadrant of the complex plane. It seems that this requirement is more restrictive. 
If f{x) is a differentiable but not analytic function, then we may construct a suitable 
Filon-type method which utilizes Hermite interpolation to compute such transforms. 
These results will be reported in our future work. 
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